%Main script for dynamicPlume
clear all; close all;
load('Topography.mat')

rhoair = 1;
rhow = 1000;

%UPPER MANTLE VISCOSITY
% mu_l = 7e19; %E = 0, constant plume visc, a = 30
% mu_l = 2e19; %E = 0, constant plume visc, lower bound
% mu_l = 2e20; %E = 0, constant plume visc, upper bound
% mu_l = 4e19; %E = 5e-19, constant plume visc, a = 30
% mu_l = 3e19; %E = 1e-18, constant plume visc, a = 30
% mu_l = 2.5e21; %E = 0, constant plume visc, a = 100
% mu_l = 1.5e23; %E = 0, constant plume visc, a = 400
% mu_l = 5e23; %E = 0, constant plume visc, a = 400, upper
% mu_l = 1e22; %E = 1e-18, constant plume visc, a = 400

% mu_l = 8e19; %rs, E = 0
% mu_l = 2e19; %rs, E = 5e-19
mu_l = 2e19; %rs, E = 1e-18


% mu_l = 1.5e20; % azores, E = 0
% mu_l = 2e19; % azores, E = 5e-19
% mu_l = 2e19; % azores, E = E-18



% PLUME VISCOSITY
% mu_u = 1e17; %E = 0, constant plume visc, a =30
% mu_u = 1e16; %E = 0, constant plume visc, lower bound
% mu_u = 1e18; %E = 0, constant plume visc, upper bound
% mu_u = 1e19; %E = 5e-19, constant plume visc, a =30
% mu_u = 3e19; %E = 1e-18, constant plume visc, a =30
% mu_u = 1e17; %E = 0, constant plume visc, a =100
% mu_u = 1e17; %E = 0, constant plume visc, a = 400
% mu_u = 1e16; %E = 0, constant plume visc, a = 400, lower
% mu_u = 1e18; %E = 0, constant plume visc, a = 400, upper
% mu_u = 1e20; %E = 1e-18, constant plume visc, a = 400


% mu_u = 1e17; %rs, E = 0
% mu_u = 1e19; %rs, E = 5e-19
mu_u = 2e19; %rs, E = 1e-18


% mu_u = 1e17; %azores, E = 0
% mu_u = 5e18; %azores, E = 0, alt
% mu_u = 1.8e19; %azores, E = 5e-19
% mu_u = 1.8e19; %azores, E = 1e-18



%PLUME FLUX km^3/Myr
% V_flux = 4e5; %E = 0, constant plume visc, a =30
% V_flux = 1e6; %E = 0, constant plume visc, upper bound
% V_flux = 1e5; %E = 0, constant plume visc, lower bound
% V_flux = 3e6; %E = 5e-19, constant plume visc, a =30
% V_flux = 6e6; %E = 1e-18, constant plume visc, a=30
% V_flux = 4e5; %E = 0, constant plume visc, a =100
% V_flux = 4e5; %E = 0, constant plume visc, for figure 6 and 4, a = 400
% V_flux = 1e6; %E = 0, constant plume visc, a = 400, lower
% V_flux = 1e5; %E = 0, constant plume visc, a = 400, upper
% V_flux = 8e6; %E = 1e-18, best fit flux, a = 400

% V_flux = 2.3e5; %rs flux E = 0
% V_flux = 2.2e6; %rs flux E = 5e-19
V_flux = 3.7e6; %rs flux E = 1e-18

% V_flux = 1.5e5; %azores E = 0
% V_flux = 1.8e6; %azores, E = 5e-19, a =30
% V_flux = 4e6; %azores E = 1e-18


% Total time plume/model is active
% t_f = 2.7e7; % iceland
t_f = 2.8e7; %rs
% t_f = 3e7; %azores

% t_f = .5e7; %ice time plots
% t_f = 1e7;
% t_f = 1.5e7;
% t_f = 2e7;
% t_f = 2.5e7;
% t_f = 3e7;

%Density of upper mantle
rho_l = 3200;

%Density of lower mantle
rho_lm = 4000;
% rho_lm = 1e20; %for R = 0 approx

%Temperature excess of plume relative to upper mantle
delTemp = 200;
% delTemp = 150; %az,new
% delTemp = 0; %no density contrast, isodense model (no bump for plume channel topo)

[T, cu, cl, r, V_u, rho_u, P_l, cr_topo] = dynamicPlume_new(mu_l, mu_u, rho_l, rho_lm, delTemp, V_flux, t_f);

%Not necessary after 1/20/23, Wiki added into main isostatic balance for
%code
%Water-load topography
% wloading = (rho_l-rhoair)/(rho_l-rhow);
% T = T.*wloading;

%%

%% Volumes

r_o = 100*1000;
r_f = 5000*1000;
delt = 1e10; %seconds
% V_flux = 1e5 * 3.154e-5; %to m^3/s
delr = 5; %kilometers
delr = delr * 1000; %km to m
a_0 = 100;
a_0 = a_0*1000;
a_0_array = a_0*ones(size(r));
speryr = 3.154e7;
t_f = speryr*2.7e7;

channel_vol_pred = 8e6 * 27; %km^3, total volume added to channel over 27 myrs

% % Act_vol = 2*pi*(((r_o)*V_o)/3.154e-4).* (t_f./speryr./1e6); % total volume flux based on initial conditions
% Act_vol = 2*pi*((r_o)*V_o).* (t_f);
% final_a = (2*pi*delr).*r.*((cu+cl)-a_0_array); % change in a layer volume (should be about act_vol)

for i  = 1:(size(cu,2)-1)
    r_mid(i) = (r(i+1)+r(i))./2;
    final_c(i) = (pi*delr).*((r(i)+delr).*cu(i+1) + r(i).*cu(i)); %m^3 volume array
end
for i  = 1:(size(cu,2)-1)
    final_a(i) = (pi*delr).*((r(i)+delr).*((cu(i+1)-cl(i+1))-a_0) + r(i).*((cu(i)-cl(i))-a_0)); %m^3 volume array
end

finalav = cumsum(final_a);
finalav_sum = sum(final_a);
finalav_sum_km = finalav_sum./(1e9);
finalav_km = finalav./(1e9);

% final_c = (pi).*(r.^2 - (r-delr).^2).*(cu); %m^3 volume array
final_c_area = (r - (r-delr)).*(cu); %m^2 area array
finalcv = cumsum(final_c); %cumulate sum of volume
finalcv_sum = sum(final_c); %cumulate sum of volume
% finalcv_sum_23 = sum(final_c(2:end)); %cumulate sum of volume
finalcv_sum_km = finalcv_sum./(1e9); %cumulate sum of volume
final_c_areasum = sum(final_c_area); %cumulate sum of volume
% final_c_areasum_23 = sum(final_c_area(3:end)); %cumulate sum of volume
finalcv_km3 = finalcv./(1e9);
finalsv = sum(final_c); %total volume of the channel

channelflux_m3_27myr = finalsv; %volume added in 27 Myr
channelflux_km3_27myr = (channelflux_m3_27myr./(1e9)); %m3 to km3
channelflux_km3_myr = channelflux_km3_27myr./27;

figure % plot volume comparison before and after of a
% plot(r./1000,finalav,'-k')
hold on
% finalcv_km3(end) = finalcv_km3(end)-((2*pi*r(1)*cu(1)/2)/1e9); %add half width volume
plot(r_mid./1000,finalcv_km3,'-g')
% plot(r./1000,finalcv_km3,'-g')
plot(r./1000,channel_vol_pred.*ones(size(r)),'-r')
% plot(r./1000,Act_vol.*ones(size(r)),'-g')  % plot constant line of anticipated volume
xlabel('Distance from Plume (km)')
ylabel('Difference in Volume of Upper Mantle Layer ''a'' (m^3)')
beautifyFigureAndMaximize();

% V_u_total  = sum(V_u);
% area_total = V_u_total*delt;


% Vf_u = V_u.*(2*pi*(r_mid(2))); %calculate volume fluxes between each grid point, volume flux at each time step between index 2 and 3
% t_arr = 0:delt:(t_f); %array of model durations at each time step
% Vol_u = Vf_u.*t_arr; %total volume added from between index 2 and 3 to end of array
% vol_total = sum(Vf_u);


%% Plotting
sub = cu*(rho_u-rho_l)./(rhow-rho_lm);

figure
hold on
plot(r./1000,T,'-k')
% plot(r./1000, (T+cr_topo), '-m')
plot(r./1000,sub,'--k')
xlim([0 2000])
xlabel('Distance from plume center (km)')
ylabel('Topography (m)')
% ylim([0 600])
% 
% plot(lenICEn,elev_s,'-b')
% plot(lenICEs,elev_reyk_s,'-r')

% plot(lenAFAR_arab,elev_arab_s,'-r')
% plot(lenAFAR_goa,elev_goa_s,'-b')
% plot(lenAFAR_nub,elev_nub_s,'-g')
% plot(lenAFAR_rs,elev_rs_s,'-m')
% plot(lenAFAR_som,elev_som_s,'-c')

plot(lenaz_2,elev_2_s,'-r')
plot(lenaz_3,elev_3_s,'-b')
plot(lenaz_5,elev_5_s,'-g')
plot(lenaz_7,elev_7_s,'-m')
plot(lenaz_9,elev_9_s,'-c')

figure
hold on
%plot(r./1000, T, '-k')

plot(r./1000, cu./1000, '-r', LineWidth=2, DisplayName='cu')
plot(r./1000, (cu+cl)./1000, '-b', LineWidth=2, DisplayName='cl+cu')
plot(r./1000, (cl)./1000, '-k', LineWidth=2, DisplayName='cl')
xlabel('Distance from plume center (km)')
ylabel('Depth below lithosphere (km)')
% axis equal
xlim([0 2000])
% ylim([0 400])
set(gca, 'YDir','reverse')
legend()

figure
hold on
semilogx(r./1000, P_l, 'k', LineWidth = 2, DisplayName='P_l')
xlabel('Distance from plume center (km)')
ylabel('P_l (Pa)')
legend()

%% Save model topography

%  save('Iceland_mu_l_high')
% save('Iceland_mu_l_low')
%  save('Iceland_mu_l_high_a400')
% save('Iceland_mu_l_low_a400')

% save('Iceland_mu_u_high')
% save('Iceland_mu_u_low')
% save('Iceland_mu_u_high_a400')
% save('Iceland_mu_u_low_a400')

% save('Iceland_v_high')
% save('Iceland_v_low')
% save('Iceland_v_high_a400')
% save('Iceland_v_low_a400')

% save('Iceland_Topo_F_0_a30')
% save('Iceland_Topo_F_5e19_a30')
% save('Iceland_Topo_F_18_a30')
% save('Iceland_Topo_F_0_a400')
% save('Iceland_Topo_F_18_a400')
% save('Iceland_Topo_F_0_a100')


% save('Afar_Topo_F_0')
% save('Afar_Topo_F_5e19')
% save('Afar_Topo_F_18')


% save('Az_Topo_F_0')
% save('Az_Topo_F_5e19')
% save('Az_Topo_F_18')



% save('Ice_F_1E19')
% save('Ice_F_5E19')
% save('Ice_F_18')
% save('Ice_F_1E19_a400')
% save('Ice_F_5E19_a400')
% save('Ice_F_18_a400')


% save('Ice_t_5')
% save('Ice_t_10')
% save('Ice_t_15')
% save('Ice_t_20')
% save('Ice_t_25')
% save('Ice_t_30')


